clear
% This script was written to study the influence of correlated
% and uncorrelated noiseon IEF and EEF on thier Coherence , Cross spectra
% and transfer function. The IEF was simulated by a power spectra with
% slope - 2 and EEF was simulated with a slope -2 spectra plus a predicted
% PP EF (using 2008 TF and simulated IEF) plus daily harmonics (6,8and12
% hours). A large number of data points used to average out the effect of
% daily harmonics in the TF.
 

% Background spectra of slope -2
%define coefficients

A = [1 -.999];
numsamples = 1024*1000;
nfft = 1024*2;
fs = 1/300;
pp_amplification = 1;

% calculate time series
background_eef = filter(1,A,randn(numsamples,1))/20;
%[Pxx_eef_bk,F_eef_bk]=pwelch(background_eef,nfft,1,nfft,(fs));
%loglog(1./(F_eef_bk*3600), abs(Pxx_eef_bk),'r');
%axis([1e-1 1e2 1e0 1e8]);
%hold on;

ief = filter(1,A,randn(numsamples,1));
ief_copy = ief;% + filter(1,A,randn(numsamples,1)); % add uncorrelated noise to show effect of bias
%[Pxx_ief_bk,F_ief_bk]=pwelch(ief,nfft,1,nfft,(fs));
%loglog(1./(F_ief_bk*3600), abs(Pxx_ief_bk),'k');



% Add synthetic PP spectra
% estimates from JULIA & CHAMP. Derived on 12/09/2008
Tx_a = [ 1.0000, -0.6023, -0.6600, 0.4451, -0.0463 ];
Tx_b = [ 0.0052, 0.0151, 0.0014, -0.0152, -0.0061 ];

pp_eef = filter(Tx_b,Tx_a, ief );
%pp_eef = filter(Tx_b,Tx_a, filter(1,A,randn(numsamples,1)) );
%pp_eef = filter(1,A,randn(numsamples,1));
%[Pxx_eef_pp,F_eef_pp]=pwelch(pp_eef,nfft,1,nfft,(fs));
%loglog(1./(F_eef_pp*3600), abs(Pxx_eef_pp),'g');

% Add daily variations

Amp = [1 0.75 0.5 0.01]/4; 
Frq = [1/24 1/12 1/8;]./3600; % 24 hur, 12 hr and 8 hour
x = (1:numsamples)/fs;
y =     Amp(1)*sin(2*pi*x*Frq(1)+pi/4) ...
    +   Amp(2)*sin(2*pi*x*Frq(2)+pi/4) ...
    +   Amp(3)*sin(2*pi*x*Frq(3)+pi/4) ...
    +   normrnd(0,Amp(4),1,length(x));
%[Pxx_eef_daily,F_eef_daily]=pwelch(y,nfft,1,nfft,(fs));
%loglog(1./(F_eef_daily*3600), abs(Pxx_eef_daily),'g');

% ALL EEF
eef = pp_eef * pp_amplification + background_eef/2 + y';

%[Pxx_eef,F_eef]=pwelch(eef,nfft,1,nfft,(fs));
%loglog(1./(F_eef*3600), abs(Pxx_eef),'r');


%ief = ief_copy;
%% calculating Coherence between IEF and EEF

%ief = ief_copy + filter(1,A,randn(numsamples,1))/2;

[Pxx,F] = pwelch(ief, hanning(nfft),1,nfft,fs);
[Pxy,F] = cpsd(ief, eef, hanning(nfft),1,nfft,fs);
[Pyy,F] = pwelch(eef, hanning(nfft),1,nfft,fs);
[Txy,F] = tfestimate(ief, eef , hanning(nfft),1,nfft,fs);
[Cxy,F] = mscohere(ief, eef, hanning(nfft),1,nfft,fs);

[Txy_exp,F] = tfestimate(ief_copy, pp_eef, hanning(nfft),1,nfft,fs);
[Cxy_exp,F] = mscohere(ief_copy, pp_eef, hanning(nfft),1,nfft,fs);


N_data = numsamples/nfft;
Err = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy)./Cxy ) ) .* abs(Txy);


%% Plot TF, Coherence

figure(1);
h1=subplot(311);
set(h1,'Position',[ 0.15 0.5 0.8 0.45]);
loglog(1./(F*3600),abs(Pxy),'k', 'LineWidth', 2);
set(gca,'FontSize',16);
hold on;
loglog(1./(F*3600),abs(Pxx),'r', 'LineWidth', 2);
loglog(1./(F*3600),abs(Pyy),'b', 'LineWidth', 2);
legend('Cross','IEF','EEF','Interpreter', 'Latex');
%xlabel('Period in Hours');
ylabel('Power Spectral Density $[\frac{mV^2}{m^2}] / Hz$', 'Interpreter', 'Latex');
axis([1e-1 1e2 1e-1 1e7]);
yticks = [1e-1 1e1 1e3 1e5];
set(gca,'YTick',yticks);
set(gca,'XTickLabel',[]);
grid on;
h2=subplot('Position',[ 0.15 0.27 0.8 0.18]);

 %loglog(1./(F*3600),abs(Txy),'r', 'LineWidth', 2);
 loglog((1./(3600*F)),abs(Txy),'ks','MarkerSize',2, 'LineWidth',1);
hold on;
 for i = 1: length(Txy);
    plot([1./(3600*F(i)) 1./(3600*F(i)) ], [abs(Txy(i)) - ... 
        Err(i) abs(Txy(i)) + Err(i)],'k','LineWidth',1);
    hold on;
end;
loglog(1./(F*3600),abs(Txy_exp)*pp_amplification,'r', 'LineWidth', 2);

 set(gca,'FontSize',16);
%legend('Transfer Function','Recov','Exp');
%xlabel('Period in Hours');
ylabel('TF Magintude','Interpreter', 'Latex');
axis([1e-1 1e2 1e-3 1e0])
yticks = [1e-3 1e-2 1e-1 1e0];
 set(gca,'YTick',yticks);
set(gca,'XTickLabel',[]);
grid on;
 set(gca,'YMinorGrid','off');
 set(gca,'XMinorGrid','off')
h3=subplot('Position',[ 0.15 0.1 0.8 0.15]);
 % semilogx(1./(F*3600), abs(Cxy_exp),'r','LineWidth', 2);
 %hold on;
 semilogx(1./(F*3600),Cxy,'k', 'LineWidth', 2);
 set(gca,'FontSize',16);
 %legend('Coherence');
xlabel('Period (Hours)','Interpreter', 'Latex');

ylabel('Coherence','Interpreter', 'Latex');
%legend('Original','Recovered');
axis([1e-1 1e2 0 1]);
yticks = [0:0.2:1];
set(gca,'YTick',yticks);
grid on;

set(gcf, 'PaperPositionMode', 'auto');
% 
%  saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/synthetic_tf_testing.fig');
%  saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/synthetic_tf_testing.png','png');
%  saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/synthetic_tf_testing.eps','eps');


%% June 21 (really ?) adding the Aliasing plot

% creating a new eef time series by down sampling the background eef by
% every 18th value (90 min interval)

eef_down = background_eef(1:18:end);

% Power spectra of the background EEF
[Pyy,F] = pwelch(background_eef, hanning(nfft),1,nfft,fs);
%Power spectra of the resampled signal
[Pyyd,Fd] = pwelch(eef_down,hanning(128),0,128,1/(90*60)); % 90 minutes interval

% now plot
loglog(1./(3600*F(2:end-1)),Pyy(2:end-1),'b','LineWidth', 2);
set(gca,'FontSize',16);
hold on;
loglog(1./(3600*Fd(2:end-1)),Pyyd(2:end-1),'r','LineWidth', 2);
ylabel('Power Spectral Density $[\frac{mV^2}{m^2}] / Hz$', 'Interpreter', 'Latex');
xlabel('Period (Hours)','Interpreter', 'Latex');

axis([1e-1 1e2 1e-1 1e5]);
yticks = [1e-1 1e1 1e3 1e5];
set(gca,'YTick',yticks);
grid on;

legend('Orignial EEF','Aliased EEF');
